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Dynamic Hubbard models are extensions of the conventional Hubbard model that take into ac- 
count the fact that atomic orbitals expand upon double occupancy. These models give rise to 
superconductivity driven by lowering of kinetic energy when the electronic energy band is almost 
full, with higher transition temperatures resulting when the ions are negatively charged. It is shown 
here that systems described by dynamic Hubbard models have a tendency to expel negative charge 
from their interior to the surface, and they have a tendency to develop charge inhomogeneity in the 
bulk. These tendencies are largest in the parameter regime where the models give rise to highest 
superconducting transition temperatures. We propose this physics as an explanation for the charge 
inhomogeneity, negatively charged grain boundaries and negatively charged vortices observed in 
high temperature superconductors. Below the superconducting transition temperature the models 
considered here describe a negatively charged superfluid and positively charged quasiparticles, un- 
like the situation in conventional BCS superconductors where quasiparticles are charge neutral on 
average. We examine the temperature dependence of the superfluid and quasiparticle charges and 
conclude that spontaneous electric fields should be observable in the interior and in the vicinity of su- 
perconducting materials described by these models at sufficiently low temperatures. We furthermore 
suggest that the dynamics of the negatively charged superffuid and positively charged quasiparticles 
in dynamic Hubbard models can provide an explanation for the Meissner effect observed in those 
superconducting materials that are not described by conventional BCS-London theory. 

PACS numbers: 



I. INTRODUCTION 

The conventional single band Hubbard model with 
Hamiltonian Jij i2j 

H = -Y^ [Ujcl^Cja + h.c] + Uj^^it^^l (1) 

i,j,<j i 

has been used to describe the physics of many real ma- 
terials. The model ignores the fact that non- degenerate 
atomic orbitals are necessarily modified by double elec- 
tronic occupancy [3^-7]. To remedy this deficiency a vari- 
ety of new Hamiltonians have been proposed and studied 
that we will generically call 'dynamic Hubbard models', 
that take into account the fact that orbital expansion 
takes place when a non-degenerate atomic orbital is dou- 
bly occupied |5H^. 

One (of several f5T]) ways to incorporate this physics is 
by the substitutionyjj i22, .23J 

Uui^Uii ->([/ + aqi)ni^nii (2) 

where a is a coupling constant (assumed positive) and 
Qi a local boson degree of freedom describing the orbital 
relaxation, with equilibrium position at = if zero 
or one electrons are present: upon double occupancy of 
the orbital at site i, qi will change from zero to a neg- 
ative value to reduce the electronic on-site repulsion, to 
an extent determined by the boson dynamics. If we de- 
scribe the dynamics of this boson by a simple harmonic 



oscillator [22] 

H^ = ^ + ^Kqf + {U + aq,)n,^n,^ (3) 

the on-site repulsion is reduced from U to Uef f = U — 
/{2K) when qi takes the value qi — —a/K correspond- 
ing to minimum energy when the site is doubly occupied, 
versus g^ = when the site is unoccupied or singly oc- 
cupied. The conventional Hubbard model does not allow 
the orbital to relax, in other words it corresponds to the 
limit if CXI of an infinitely stiff orbital. 

The Hamiltonian Eq. (3) breaks electron-hole 
symmetry [TT]: the importance of this physics increases 
as the filling of the electronic energy band increases, 
as is simply seen by taking the mean field expectation 
value of Eq. (3). In addition, the importance of this 
physics increases when the ionic charge is small^llj. since 
in that case the orbital expansion is larger (for exam- 
ple, the orbital expansion is larger for H~ than for He), 
corresponding to a smaller stiffness parameter K in Eq. 
(3). These two facts imply that the importance of this 
physics increases the more negative charge the system 
has. Thus, it will come as no surprise that a system 
where this physics is important will have a strong ten- 
dency to expel negative charae ^I^ I25j . This expulsion 
occurs already at the atomic level of course since the or- 
bital expansion upon double occupancy implies radially 
outward motion of negative charge. 

When an atomic orbital expands, the electronic kinetic 
energy is lowered, since in an orbital of radial extent r 
the electron kinetic energy is of order / (2mi,r'^) , with 
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nie the electron mass. Thus, a system described by a dy- 
namic Hubbard model will have lower electronic kinetic 
energy than if the orbital relaxation effect didn't happen. 
Just like at the atomic level the orbital expansion and 
associated negative charge expulsion is associated with 
lowering of kinetic energy, we will find for the system 
as a whole described by a dynamic Hubbard model that 
negative charge expulsion is associated with lowering of 
kinetic energy. Note however that our site Hamiltonian 
Eq. (3), describing orbital expansion when the value of 
Qi is negative, does not have a term explicitly describing 
the kinetic energy lowering of the atomic electron when 
Qi adopts a negative value. 

There have been theoretical suggestions in the litera- 
ture that electronic kinetic energy lowering may provide 
an unconventional mechanism for suDerconductivitv |10l 
[26U43j . applicable to materials that are not described 
by the conventional BCS-London theory that predicts 
an increase in kinetic energy upon the transition to 
superconductivity |44| . There have also been experi- 
mental indications that in certain materials and doping 
regimes the transition to superconductivity may be as- 
sociated with kinetic energy lowering |45H54j . Since dy- 
namic Hubbard models describe kinetic energy lowering 
a natural question is whether these models can describe 
superconductivity in cases where superconductivity is as- 
sociated with lowering of kinetic energy. This is indeed 
the case [55H58] . 

In addition, dynamic Hubbard models are by nature 
electron-hole asymmetric and so are superconductors, as 
evidenced by the fact that a rotating superconductor de- 
velops a magnetic field that is always parallel, never an- 
tiparallel, to the direction of the mechanical angular mo- 
mentum of the bodv[5^. This suggests that dynamic 
Hubbard models may be more appropriate to describe 
superconductors than the conventional Hubbard model 
that is electron-hole symmetric. 

As a consequence of the charge expulsion physics the 
superconducting state in systems described by dynamic 
Hubbard models has quasiparticles that are positively 
charged on average [60], and as a consequence the super- 
fluid has excess negative charge, in contrast to conven- 
tional BCS-London superconductors where quasiparticles 
are charge neutral on average. We will examine the con- 
sequences of this physics for superconductors described 
by these models at temperatures well below the super- 
conducting transition temperature. 

The potential (Coulomb) energy of a system of charges 
is minimized when the charge distribution is uniform. A 
non-uniform charge distribution gives rise to electrostatic 
fields and an associated potential energy cost. It will be 
favored if this cost is compensated by a kinetic energy 
gain, i.e. lowering of kinetic energy. Thus it is natural 
to expect, as we will find in this paper, that dynamic 
Hubbard models are prone to develop charge inhomo- 
geneity, and in extreme cases charge separation, where 
the kinetic energy lowering overcompensates for the po- 
tential energy cost. Many materials of recent interest, in- 



cluding high Tc cuprates and manganites, exhibit charge 
inhomogeneitv pTM65j . suggesting that dynamic Hubbard 
models may be useful to describe them. 

In superconductors described by the conventional 
BCS-London theory, no negative charge expulsion oc- 
curs, the kinetic energy is raised rather than lowered 
in the transition to superconductivity, quasiparticles are 
charge neutral on average, and the Meissner effect is 
argued to be completely understood within the frame- 
work of the conventional theory [55H76j . In this paper we 
do not attempt to describe superconductors described 
by the conventional theory by dynamic Hubbard mod- 
els. Rather, our aim is to describe superconductors that 
are not described by conventional BCS-London theory. 
There is general agreement among workers in the field 
that there are several classes of superconducting materi- 
als that are not described by conventional theorv [77ll78j . 
All those materials exhibit a Meissner effect, and as men- 
tioned there is experimental evidence for at least some 
unconventional materials that the kinetic energy is low- 
ered in the transition to superconductivity. Since those 
materials are not described by the conventional theory 
it is certainly not obvious that the explanation of the 
Meissner effect of conventional superconductors provided 
by the conventional theorv [5SH75] would apply to these 
non-conventional materials. The negative charge expul- 
sion physics driven by kinetic energy lowering of dy- 
namic Hubbard models provides a natural explanation 
for the Meissner effect observed in these unconventional 
materials |79i i8U]: just like in classical plasmas obeying 
Alfven's theorem [ST], the magnetic field lines move with 
the expelled negative charge. 

II. DYNAMIC HUBBARD MODELS 

Figure 1 depicts the essential physics of dynamic Hub- 
bard models as opposed to conventional Hubbard mod- 
els: when a second electron (of opposite spin) is added 
to a non-degenerate orbital, it expands. We can describe 
this physics by a multi-orbital tight binding model (at 
least two orbitals per site) [HI [H2 , or with a background 
spin[31 [TU] or harmonic oscillator [S31 [Hlj degree of free- 
dom that is coupled to the electronic double occupancy, 
as in Eq. (2). Assuming the latter, the site Hamiltonian 
is given by Eq. (3), and the Hamiltonian can be written 
as 

iJ = - ^ [ttj4„Cja + h.c] + ^ huoala^ 

+ ^[U + ghLUo){al + ai)]ni^nii (4) 

i 

with frequency cjq = K/m and g — a/{2KfkjJoY^'^ the 
dimensionless coupling constant. Estimates for the val- 
ues of these parameters were discussed in ref. [11]. In 
particular, for Is orbitals g'^uJo ^ A.leV . Quite generally 
we expect g to increase proportionally to 1 / Z and ujq to 
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FIG. 1: In the conventional Hubbard model the atomic or- 
bital is not modified by electronic occupancy. In the dynamic 
Hubbard model and in real atoms, addition of the second 
electron causes orbital expansion due to the electron-electron 
interaction. Negative charge is expelled outward and the ki- 
netic energy of the electrons is lowered relative to that with 
a non-expanded orbital. 



increase proportionally to Z^, where Z is the charge of 
the ion when the orbital under consideration is empty [TT]. 
However, in an even more realistic model ujq should also 
change with different electronic occupation. That issue 
is beyond the scope of this paper. 

Using a generalized Lang-Firsov transformation [Tl] 
[55H57] the electron creation operator cj^ is written in 
terms of new quasiparticle operators c] as 



J — p5(al-ai)ni,-^st 
So- ~ Ha 



+ iii^-a X {incoherent part) (5) 



where the incoherent part describes the processes where 
the boson goes into an excited state when the electron is 
created at the site. For large wq those terms become small 
and can be neglected, and even for not so large loq we 
have found from numerical studies that their effect does 
not change the low energy physics qualitativelv|881 155]. 
Hence we will ignore those terms in what follows, which 
amounts to keeping only ground state to ground state 
transitions of the boson field. 

The electron creation operator is then given by 



(6a) 



so that it decreases monotonically from 1 when the band 
is almost empty to 5^ < 1 when the band is almost full. 
The single particle Green's function and associated spec- 
tral function is renormalized by the multiplicative factors 
on the quasiparticle operators given in Eq. f 6a)) .8(1), 
which on the average amounts to multiplication of the 
spectral function by the quasiparticle weight Eq. (6c). 
This will cause a reduction in the photoemission spectral 
weight at low energies from what would naively follow 
from the low energy effective Hamiltonian, an effect ex- 
tensively discussed in Ref. [5S| and recently rediscovered 
in 20 . A corresponding reduction occurs in the two- 
particle Green's function and associated low frequency 
optical properties [86l l88|. 

The low energy effective Hamiltonian is then 



cLcja + h.c] + Ueff fli^fltl (7a) 



- [1 + (5 - l)n.,-.][l + iS- l)n„ _,]t. 



(7b) 



and Ueff = U — huj^g^. Thus, the hopping amplitude 
for an electron between sites i and j is given by tij , Stij 
and S'^tij depending on whether there are 0, 1 or 2 other 
electrons of opposite spin at the two sites involved in the 
hopping process. 

The physics of these models is determined by the mag- 
nitude of the parameter S, which can be understood as 
the overlap matrix element between the expanded and 
unexpanded orbital in Fig. 1. It depends crucially on the 
net ionic charge Z, defined as the ionic charge when the 
orbital in question is unoccupied In Fig. 1, Z = 1 ii 
the states depicted correspond to the hydrogen ions , 
H and and Z = 2 if they correspond to He~^^, ffe"*" 
and He. In a lattice of anions, as in the Cu — O 
planes of high Tc cuprates, the states under considera- 
tion are O, 0~ and and Z = 0, and in the B~ planes 
of MgB2, Z = 1. The effects under consideration here 
become larger when S is small, hence when Z is small. 
An approximate calculation of S as function of Z is given 
in[II]. 

We now perform a particle-hole transformation since 
we will be interested in the regime of low hole concentra- 
tion. The hole creation operator is given by, instead of 
Eq. (6a) 



(8a) 



where fii^a is now the hole site occupation, and the hole 
quasiparticle weight increases with hole occupation n as 



S = e-s'/2 



(6b) 



and the quasiparticle weight for electronic band filling n 
(n electrons per site) is 



z(n) = il + iS-l)-r 



3c) 



z,(n) = 5^(1 + (1-1)^)2 



(8b) 



For simplicity of notation we denote the hole creation 
operators again by c|^, the hole site occupation by ni^ 
and the effective on-site repulsion between holes of oppo- 
site spin Ueff (the same as between electrons) by U to 
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simplify the notation. The Hamiltonian for holes is then 

+ h.c] + uY, ^^T'^a (9a) 

ija i 



with t'^j = S^tij the hopping amplitude for a single hole 
when there are no other holes in the two sites involved 
in the hopping process. The hole hopping amplitude 
and the effective bandwidth increase as the hole occupa- 
tion increases, and so does the quasiparticle (quasihole) 
weight Eq. (8b). 

Finally, we will assume there is only nearest neigh- 
bor hopping tij = t for simplicity and write the nearest 
neighbor hopping amplitude resulting from Eq. (9b) as 

tfj = </i + M(ni-a + f^J,-c^) + ^t2ni^_anj-a (10a) 

with 

th = tS^ (10b) 
M = tS{l - S) (10c) 

Ai2 = t{l - Sf = {Atf/th. (lOd) 

The non-linear term with coefficient At2 is expected to 
have a small effect when the band is close to full (with 
electrons) and is often neglected. Without that term, 
the model is also called generalized Hubbard model or 
Hubbard model with correlated hopping [501 [OT] . The ef- 
fective hopping amplitude for average site occupation n 
is, from Eq. (lOa) 

2 

Tl 

t{n) =th + nAt+ —At2 (11) 

so that a key consequence of integrating out the higher 
energy degrees of freedom is to renormalize the hopping 
amplitude and hence the bandwidth and the effective 
mass (inverse of hopping amplitude), as emphasized since 
the beginning of our studies of these Hamiltonians [H2 
and recently rediscovered in a different context [20] . 

III. GENERALIZED DYNAMIC HUBBARD 
MODELS 

More generally, one can assume that the boson degree 
of freedom will couple not only to the double orbital oc- 
cupancy but also to the singly occupied orbital [55]. The 
site Hamiltonian is then 

Hi = fiwoa\ai + hwolgui^n^i + go{n^-^ + n^)] (12) 



At the atomic level, the coupling go will appear when con- 
sidering an orbital for atoms other than hydrogenic ones, 
and represents the modification of the states of electrons 
in other orbitals in the atom when the first electron is 
created in the orbital under consideration. We expect 
this effect to be much smaller than the modification of 
the state of the electron residing in the same orbital, the 
physics described by g, hence go << g, particularly when 
the ionic charge Z is small. The formal development for 
the site Hamiltonian given by Eq. (12) is very similar to 
the one discussed in the previous section and is given in 
Ref. [SB]. In particular, Eq. (6) becomes 

cL = [r+(5-T)n,_^]4 (13a) 
T = (13b) 

z(n) = (r+(5-T)^)^ (13d) 

and of course S < T always [10, since g > 0. The hop- 
ping amplitudes given in the previous section are simi- 
larly modified by replacing 1 by T in various places, as 
discussed in Ref. |86| . go gives a renormalization of the 
quasiparticle mass, bandwidth and quasiparticle weight 
that is independent of band filling, and g gives a band- 
filling dependent contribution. 

Recently, a single band model with site Hamiltonian 
of the form Eq. (12) with g = and go was 
considered [20] to describe the effect of higher energy elec- 
tronic excitations on the low energy electronic physics 
within dynamical mean field theory]^. In our view 
this is an unphysical limit since we expect g >> go 
quite generally. Some of the effects discussed in refs. 
[13 [53 IHSl [Ml IS2 were rediscovered in that work[10], in 
particular the "bandwidth renormalization from Coulom- 
bic screening" , the "renormalization of the hopping and 
hybridization parameters", and the fact that "the pho- 
toemission spectral weight in the frequency range de- 
scribed by the effective model is reduced" , as well as the 
fact that the effective model "gives a reasonable good 
description" even if wo is not too large[2D]. The effects 
specifically associated with g will presumably be re- 
discovered when the authors of Ref. [2D] include g in their 
Hamiltonian, as well as the effect on two-particle Green's 
functions and optical properties discussed in [TD], [2T], 
[52, [HI] and [55], when the authors of Ref. [10] consider 
those quantities. 

IV. HOLE PAIRING AND 
SUPERCONDUCTIVITY IN DYNAMIC 
HUBBARD MODELS 

As we[5S] and others [51H57] have discussed, the corre- 
lated hopping At gives a strong tendency to pairing and 
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FIG. 2: Expected charge distribution in the ground 
state of superconductors according to the theory of hole 
superconductivity |25] . 



superconductivity when a band is almost full. The hop- 
ping amplitude for a single hole is tf^, and it increases to 
tfi + At when the hole hops to or from a site occupied 
by another hole (of opposite spin), thus giving an incen- 
tive for holes to pair to lower their kinetic energy. The 
superconductivity described by this model has a num- 
ber of interesting features [551 that we have proposed are 
relevant to the description of high Tc cuprates, namely 
strong dependence of Tc on hole concentration[55] . energy 
dependent gap function and resulting tunneling asymme- 
try of universal sign|98]. superconductivity driven by ki- 
netic energy lowering and associated low energy optical 
sum rule violation|56j, change in optical spectral weight 
at frequencies much higher than the superconducting en- 
ergy gap upon onset of superconductivity[10] , strong pos- 
itive pressure dependence of r^i55j. increased quasiparti- 
cle weight upon entering the superconducting state [5S]. 
etc. Many of these predictions are supported by obser- 
vations on high Tc cuprates made both before and after 
the predictions were made. 

The more fundamental dynamic Hubbard model from 
which the At interaction derives has also been stud- 
ied using Eliashberg theorylSH] and exact numerical 
methods [13 and shows an even stronger tendency to 
pairing and superconductivity. 



V. NEGATIVE CHARGE EXPULSION IN 
DYNAMIC HUBBARD MODELS 

Within the theory of hole superconductivity, supercon- 
ductors are predicted to have a non-homogeneous charge 
distribution in the ground state [24], as shown schemat- 
ically in Fig. 2: excess positive charge in the interior 
and excess negative charge near the surface, resembling 
a "giant atom" [99] . This charge distribution is predicted 
from modified London-like electrodynamic equations [25j. 
As we will see here, dynamic Hubbard models predict 
qualitatively similar behavior. 



We consider the Hamiltonian for holes Eq. (9), with 
the hopping amplitudes given by Eq. (10), which we 
reproduce here for convenience: 

H = -Y1 + h.c] + UY, "^t"4 (14a) 

ij(7 i 

tij =th + At(rii _o. -I- nj,^cr) + At2ni^^crnj-a (14b) 

It is clear from the form of this Hamiltonian that the ki- 
netic energy decreases when the number of holes in the 
band increases, since the hopping amphtudes Eq. (14b) 
increase with hole occupation. This suggest that the sys- 
tem will have a tendency to expel electrons from its in- 
terior to the surface, because the coordination of sites in 
the interior is larger than of sites at the surface. In what 
follows we study this physics numerically. 

We assume a cylindrical geometry of radius R and in- 
finite length in the z direction. We decouple the inter- 
action terms within a simple mean field approximation 
assuming < >— ni/2 with ni the hole occupation at 
site i, and obtain the mean field Hamiltonian 

(15a) 



<ij>.,a 

(15b) 

Hmf,pot = J^n^ (15c) 

i 

Assuming a band filling of n holes per site, we diagonal- 
ize the Hamiltonian Eq. (15) with initial values Ui = n 
and fill the lowest energy levels until the occupation n is 
achieved. From the Slater determinant of that state we 
obtain new values of rii for each site, and iterate this pro- 
cedure until self-consistently is achieved. We can extend 
this procedure to finite temperatures, iterating to self- 
consistency for a given chemical potential fi. We consider 
then the resulting occupation of the sites as function of 
the distance r to the center of the cylinder. Sometimes 
there are non-equivalent sites at the same distance from 
the origin (e.g. (5,0) and (3,4)) that yield somewhat dif- 
ferent occupation, for those cases we show the average 
and standard deviation as error bars in the graphs. 

Figure 3 shows a typical example of the behavior found. 
Here we assumed At2 = 0, corresponding to the sim- 
pler Hubbard model with correlated hopping and no six- 
fermion operator term. Even for At = the hole occu- 
pation is somewhat larger in the interior than near the 
surface. When the interaction At is turned on, the hole 
occupation increases in the interior and decreases near 
the surface. This indicates that the system expels elec- 
trons from the interior to the surface. The effect becomes 
more pronounced when At is increased. 
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FIG. 3: Hole site occupation per spin for a cylinder of radius 
7i = 11 as function of r/R, with r the distance to the center, 
for a cubic lattice of side length 1. There are 377 sites in a 
cross-sectional area {nR^ = 380.1). The average occupation 
(both spins) is n = 0.126 holes per site and the temperature 
is feflT = 0.02. 
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FIG. 4: The diameters of the circles are proportional to the 
hole occupation of the site. Note that for finite At the hole 
occupation increases in the interior and is depleted near the 
surface. The parameters correspond to the cases shown in 
Fig. 3. 



FIG. 5: Kinetic, potential and total energy per site for 
At = 0.25 as function of number of iterations starting with a 
uniform hole distribution. 
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FIG. 6: Effect of Coulomb repulsion, th = 0.1, At = 0.25. As 
the on-site repulsion increases, the charge expulsion decreases 
and the occupation becomes more uniform. 



Figure 4 shows the hole site occupations as circles of 
diameter proportional to it, for the cases At = and 
At = 0.25 of Fig. 3. Note that the interior hole occupa- 
tion is larger for At = 0.25 than it is for At = 0, while 
near the surface the hole occupation is larger for At = 0. 
Again this shows that the system with At — 0.25 is ex- 
pelling electrons from the interior to the surface, thus 
depleting the hole occupation near the surface. 

These results are obtained by iteration. Fig. 5 shows 
the behavior of the energies as a function of iteration 
number for the cases At = and At = 0.25 of Fig. 
3. The initial values correspond to a uniform hole dis- 
tribution with each site having the average occupation. 
The evolution is non-monotonic because in the interme- 
diate steps the overall hole concentration increases, nev- 
ertheless it can be seen that for the case At = 0.25 the 
final kinetic energy when self-consistency is achieved is 



lower, and the final potential energy is higher, associ- 
ated with the larger hole concentration in the interior 
and the lower hole concentration near the surface shown 
in Fig. 4. This is of course what is expected. For the 
case At = instead there is essentially no difference in 
the energies between the initial uniform state and the 
final self-consistent state. 

The charge expulsion caused by At is counteracted by 
the effect of Coulomb repulsion, just like the tendency to 
pairing is [55]. Figure 6 shows the effect of increasing the 
on-site repulsion for a fixed value of At. 

The effect of the non-linear occupation-dependent hop- 
ping term At2 (Eq. (lOd)) is shown in Fig. 7. It 
increases the charge expulsion tendency relative to the 
model where this term is omitted (Hubbard model with 
correlated hopping). 

As the correlated hopping amplitude At increases, and 
even more so in the presence of At2, the system appears 
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FIG. 7: Effect of non-linear occupation-dependent hopping 
term At2, Eq. (7c) for the case At = 0.25 of Fig. 3. 
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FIG. 8: As the correlated hopping terms increase, the system 
develops a tendency to phase separation, where essentially all 
the holes condense to the interior. Parameters are the same as 
in Fig. 3 except as indicated. The maximum hole occupation 
per spin is 0.128 and 0.214 for the left and right panel, the 
average hole occupation per spin is 0.063. 




At=0 .3 , At2= At=0 .25 , At2?^ 



FIG. 9: Hole site occupations for the same parameters as Fig. 
8 but average hole occupation per spin 0.175. The maximum 
hole occupation per spin is 0.107 for the left panel, 0.116 for 
the right panel. 

corresponds to smaller S, the overlap of the atomic or- 
bitals when one and two electrons are at the orbital. As 
discussed earlier, S is smaller when the ionic charge Z is 
smaller, corresponding to a more negatively charged ion. 
The fact that the effective Hamiltonian derived from this 
physics expels more negative charge the more negatively 
charged the ion is makes of course a lot of sense and 
can be regarded as an internal consistency check on the 
validity of the model. 

For a normal metal, the charge expulsion physics will 
be compensated to a large extent by other effects since no 
electric field can exist in the interior of a normal metal. 
The Coulomb repulsion will do this, but as we argue 
in the next sections some residual effects of charge ex- 
pulsion can be seen even in the normal state. For the 
superconducting state, we have proposed new electrody- 
namic equations that give rise to "charge rigidity" |100| 
and the inability of the superfluid to screen interior elec- 
tric fields so that the charge expulsion physics can man- 
ifest itself [IS] . 



to develop a tendency to phase separation, where holes 
condense in the interior and the outer region of the cylin- 
der becomes essentially empty of holes. This happens 
very rapidly as function of the parameters for the finite 
system under consideration. Examples are shown in Fig. 
8. We return to this point in the next section. 

The charge expulsion tendency and associated effects 
described in this and other sections become rapidly 
weaker when the hole concentration increases. To illus- 
trate this we show in Fig. 9 the charge distribution for 
the same Hamiltonian parameters as Fig. 8 but average 
hole occupation per site n = 0.35 instead of n = 0.126. 
It can be seen that the charge expulsion is much smaller. 

In summary, we have seen in this section that the dy- 
namic Hubbard model promotes expulsion of negative 
charge from the interior to the surface of the system when 
the band is almost full, and that the charge expulsion 
physics is associated with kinetic energy lowering, just as 
in the single atom, Fig. 1. The charge expulsion tendency 
is largest when the parameter At is largest, which in turn 



VI. PHASE SEPARATION 

The tendency to charge expulsion in the dynamic Hub- 
bard model can also lead to a thermodynamic instability 
and ensuing phase separation. The question of phase sep- 
aration in these models in certain parameter ranges was 
recently considered by Montorsi and coworkers 1011 [TU^ . 

Let us consider the correlated hopping model first for 
simplicity {At2 — 0). The effective hopping amplitude 
for a hole is t{n) = th + nAt with n the hole density 
per site. In the dilute limit, the kinetic energy of a hole 
is e = —zt{n), with z the number of nearest neighbors 
to a site, and it is clear that when At is much larger 
than til putting for example all holes into half the sys- 
tem and thus emptying the other half will double n and 
hence decrease the kinetic energy per hole as well as the 
total energy if it is dominated by kinetic energy. This 
tendency to phase separation will be countered both by 
Pauli exclusion and by Coulomb repulsion. 
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Consider a flat density of states for simplicity. The 
effective bandwidth increases hnearly with n 



D{n) = Dh + Kn 



(16) 



with D = 2zt, Dh = 2zth, K = 2zM. The density of 
states per site per spin is given by 5(e) = \/D and the 
ground state kinetic energy by 



eg{e)de^^{n^-2n) (17) 



D/2 

with ji = {n ~ l)D{n)/2 the chemical potential for n 
holes per site. Adding the on-site repulsion in a mean 
field approximation yields 



Eo{n) 



Dh +nK 



{n^ - 2n) + 



V 



(18) 



and the system will be unstable towards phase separation 
into hole-rich and hole-poor regions when the condition 



is satisfied, hence 



or equivalcntly 



2 



X(-n-l)<0 



K > 



U + Dh 
2(1 -fn) 



At > 



th + U/{2z) 
2(1 



in) 



(19) 



(20) 



(21) 



For the parameters used in the previous section, e.g. 
tfi = 0.1, n = 0.126, U = 2 and z = 4 appropriate to 
two dimensions Eq. (21) yields At > 0.216. The ten- 
dency to phase separation becomes even stronger when 
the nonlinear term At2 is included. After some simple 
algebra Eq. (21) is modified to the condition 



At > 



4 + U/i2z) 3Ai2n(l - n)/{4th) 



2(1 



2(1 



(22) 



which for the parameters given above yields At > 0.182. 
These estimates are consistent with the numerical results 
shown in the previous section. Note that as n increases 
larger At is needed for phase separation. 

Note that the instability criterion Eq. (19) appears to 
be different from the usual criterion 



(9/i 
dn 



< 



(23) 



if fj, is given by the expression given right after Eq. (17). 
The reason is that the fi in Eq. (17) is not the true 
chemical potential but an effective one. The true chemi- 
cal potential is modified by contributions from both the 
Hubbard repulsion and the density dependent hopping 
terms. For the case with At2 = it is given by 



M = Me// + 



Un 



K 



< c]cja + h.C. > 



(24a) 



^^eff - 



Dh + nK 



(n-1) 



< c, 



. lb 

+ h.C. n(l - -) 



(24b) 



(24c) 



where the expectation value for the bond charge Eq. 
(24c) follows from Eq. (17). Hence we obtain from Eq. 

(24) 



^ u + Dh 

dn 2 



K{l--n) 



(25) 



in agreement with Eq. (19). The instability criteria Eq. 
(23) or Eq. (19) with the free energy replacing Eq, are 
also valid at finite temperature of course. 

In a real material in the normal state phase separa- 
tion will not occur because of the effect of longer-range 
Coulomb interactions. However this physics will clearly 
favor charge inhomogeneity, i.e. hole-rich regions that 
benefit from the lowering of kinetic energy acquired by 
increasing the hole concentration and thereby broadening 
the (local) bandwidth, and hole-poor regions where the 
kinetic energy cost is mitigated by narrowing of the local 
bandwidth. The shape of these regions will depend on the 
particular details of the system and merits further inves- 
tigation. A possible geometry for the hole-rich and hole- 
poor regions could be one-dimensional, i.e. stripes '62J. 
Other geometries like patches are also possible [64, 65]. 
Such charge inhomogeneities are commonly seen in high 
Tc superconductors [11] where the physics discussed here 
should be dominant. 



VII. CHARGE INHOMOGENEITY 

High Tc cuprates show a high tendency to charge 
inhomogeneity [nil [55M5] . We suggest that a dynamic 
Hubbard model may be relevant to describe this physics: 
because kinetic energy dominates the physics of the dy- 
namic Hubbard model, the system will develop charge in- 
homogeneity at a cost in potential energy if it can thereby 
lower its kinetic energy more, unlike models where the 
dominant physics is potential (correlation) energy like 
the conventional Hubbard model. 

Let us consider the effect of local potential varia- 
tions. We have shown before that in the superconduct- 
ing state of the model with At there is great sensitivity 
to local potential variations due to the slope of the gap 
function |103| . Here we find that the model is also sensi- 
tive to local potential variations in the normal state. 

We assume there are impurities in the system that 
change the local potential at some sites, and compare the 
effect of such perturbations for the dynamic and conven- 
tional Hubbard models. As an example we take parame- 
ters tfi — 0.1, U — 2 and consider site impurity potentials 
of magnitude ±0.2 at several sites as indicated in the cap- 
tion of Fig. 9. For the dynamic Hubbard model we take 
At = 0.2, At2 = 0.4, corresponding to 5 = 0.333. 
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FIG. 10: Hole site occupation per spin in a system of radius 
i? = 11 with 5 impurities at positions (-1,0), (2,2), (3,-4), (- 
5, -5), (-6, 7) with potential strength -0.2, -1-0.2, -0.2, +0.2, 
-0.2 respectively. Note the much larger variation in densities 
generated in the dynamic Hubbard model (lower panel, At2 7^ 
0) than in the conventional Hubbard model (upper panel). 
Average hole occupation per site is n = 0.126. 



Figure 10 shows the effect of these impurities on the 
charge occupation for the conventional and dynamic 
models. In the conventional Hubbard model the occu- 
pation changes at the site of the impurity potential and 
only very slightly at neighboring sites. In the dynamic 
Hubbard model the local occupation change at the im- 
purity site itself is much larger than in the conventional 
model, and in addition, the occupations change at many 
other sites in the vicinity of the impurities, as seen in the 
lower panel of Figure 10. Figure 11 shows the real space 
distribution of these changes. 

Similarly we can consider impurities where the atomic 
value of the deformation parameter S is different than in 
the bulk. This would occur for example by substituting 
an ion by another ion with different ionic charge, hence 
different orbital rigidity. For example, substituting O 
by F~ would make the orbital more rigid and increase 
S at this site, while substituting O by N would 
make the orbital more floppy and decrease S. Figure 12 
shows 5 impurities at the same locations as in Figure 9, 
with values S = 0.5 and 0.2 at the impurity sites instead 
of the bulk value S = 1/3. The larger (smaller) S will in- 
crease (decrease) the occupation. Compared to the case 
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FIG. 11: The site occupations for the case of Fig. 10, with the 
diameters of the circles proportional to the hole occupation 
of the sites. Note the 5 impurity sites at positions listed in 
the caption of Fig. 9 (three with negative potential, hence 
higher hole concentration) and two with positive potential, 
hence lower hole concentration. Note that for At — only 
the occupation at the impurity site changes appreciably, while 
for At 7^ an impurity potential of the same strength causes 
a much larger change of the occupation at the impurity site 
and occupation change also at the nearest and next nearest 
neighbor sites. 
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FIG. 12: Hole site occupation per spin in a system of radius 
i? = 11 with 5 impurities at positions (-1,0), (2,2), (3,-4), (-5, 
-5), (-6, 7) with S factor 0.5, 0.2, 0.5, 0.2, 0.5 respectively. 
All other sites have S = 0.333. n = 0.126. Note that the 
variations in density occur for even more sites than when the 
local potential is changed, Fig. 10 lower panel, even though 
the change in occupation at the impurity site itself is smaller. 



of Fig. 10, it can be seen that the change in occupation at 
the impurity site itself is somewhat smaller for these pa- 
rameters but the changes are larger at neighboring sites 
and extend to sites farther away. Similarly as in Fig. 11 
we show the real space changes in Figure 13, compared 
to the conventional Hubbard model where no change at 
all would occur since it is not sensitive to the rigidity of 
the orbital. 

The changes in site energies in the example of Figs. 10 
and 11 could occur both if there are substitutional impu- 
rities in the conducting plane with different on-site ener- 



10 



JOOOOOOOOGt 

■ ODODDOOQOOOt 

JOOOOOOQOOOOO t 
■ OOOOOOOOOOOOOOOOO 1 
iOOOOOOOOOOOOOOOOOc 

I ooooooooooooooooooo- 

jOOOOOOOOOOOOOOOOOOOt 
'OOOOOOOOOOOOOOOOOOOOC 

lODOOOOOOOOOOOOOOOQOc 
DOOOOOOOOOODOOOOOOOaO 

jooooaoooooooooooooot 

jQOaoOOOOOOQOOOOOOOOt 
OOOOOOOOOOOOOt 
' OOOOOOOOOOOOOOOOOOO 1 
lOOOOOOOOOOOOOOOOOc 
' OOOOOOOOOOOOOOOOO t 
' OOOOOOOOOOOOOOO t 
'OCOOOOOOOOOOO- 
JOOOOOOOOOt 




At=0 



At=0.2, AtjifO 



FIG. 13: The site occupations for the case of Fig. 11, with 
the diameters of the circles proportional to the hole occupa- 
tion of the sites. Same 5 impurity sites at positions listed 
in the caption of Fig. 10 (three with larger S, hence higher 
hole concentration) and two with smaller S, hence lower hole 
concentration. Note that the range of sites where the occu- 
pation changes is even larger than in Fig. 11, extending to 
third nearest neighbors. The left panel shows for comparison 
a uniform system, corresponding to a conventional Hubbard 
model that does not take into account the deformation of the 
orbital, hence S = 1 for all sites. 



gies, or if there are impurities off the plane that change 
the local potential in the plane. In contrast, the example 
of Figs. 12 and 13 would be appropriate to describe only 
impurities in the plane itself where the electrons conduct. 

The reason for this large sensitivity to local perturba- 
tions can be understood from the form of the hopping 
amplitude Eq. (14b). A change in the local occupa- 
tion will also change the hopping amplitude of a hole 
between that site and neighboring sites, which in turn 
will change the occupation of neighboring sites, and so 
on. Similarly a change in the deformation parameter S 
at a site will affect the hopping amplitudes between that 
site and its nearest neighbors, hence the occupation of 
the site and its neighbors, etc. In that way a local per- 
turbation in the dynamic Hubbard model gets amplified 
and expanded to its neighboring region, and it is easy 
to understand how a non-perfect crystal will easily de- 
velop patches of charge inhomogeneity in the presence 
of small perturbations. These inhomogeneities cost po- 
tential (electrostatic) energy, but are advantageous in ki- 
netic energy. The conventional Hubbard model does not 
exhibit this physics. 



VIII. GRAIN BOUNDARIES IN HIGH n 
CUPRATES 

The high superconducting transition temperature of 
the cuprates indicate that, if they are described by dy- 
namic Hubbard models, they have a small value of the 
parameter S or equivalently a large value of the corre- 
lated hopping parameter At. This is to be expected be- 
cause of the small ionic charge associated with the neg- 
ative ion . According to the results of the previous 



sections, this should lead to a strong tendency of high 
Tc cuprates to expel negative charge from their interior 
and to develop charge inhomogeneity in the presence of 
disorder. In this section we argue that these features 
may also be relevant to the understanding of properties 
of grain boundaries in high Tc cuprates [104H106j . 

Babcock et al |107) report results of EELS experiments 
indicating severe hole depletion in YBCO near grain 
boundaries for large angle grain boundaries, with the 
hole deficient region extending up to lOOA or more into 
the crystal. In contrast, small angle grain boundaries 
show significantly less hole depletion. It is reasonable 
to assume that for larger angle grain boundaries there 
is weaker coupling betwen the grains, and this is con- 
firmed experimentally by measurement of the Josephson 
critical current across the grain boundary |107j. Further- 
more, Babcock et al found that the structural changes as- 
sociated with the grain boundaries (structural perturba- 
tions and cation nonstochiometry) extended only about 
5A from the grain boundary into the crystal, and that the 
hole depletion is not associated with particular specimen 
preparation procedures such as time duration of oxygen 
annealing and storage conditions. The fact that the hole 
depletion region extends over a much wider region that 
associated with structural changes suggests an intrinsic 
purely electronic origin for the effect. Other workers have 
found similar results, including Browning et al |108j and 
Schneider et al[109]. 

Furthermore and consistent with this picture, it has 
been found that substituting y by Ca near grain bound- 
aries is an efficient way to prevent the hole depletion 
phenomenon'110 H112j . since Ca~^^ ions will donate fewer 
electrons to the CuO planes than the ions they 

substitute, thus increasing the hole concentration and 
thereby improving the transport properties across the 
grain boundary. 

Previous theoretical explanations of these effects have 
implicitly or explicitly assumed that grain boundaries in 
high Tc cuprates have an intrinsic positive charge that 
leads to band bending and consequently a flow of con- 
duction electrons to the vicinity of the grain boundaries 
that causes hole depletion |1 1 3H1 1 5] . However, these ex- 
planations are directly contradicted by experiments that 
measure the electric potential at the grain boundary by 
electron-beam holography [120]. These experiments show 
unequivocally that the electrostatic potential at the grain 
boundary is negative with respect to the interior [1171 - 

ns]. 

To make sense of the electron holography results, 
Mannhart suggested ^lOG that the negative potential at 
the grain boundary core may cause overdoping of the 
lower Hubbard band with holes, resulting in an empty 
band and insulating behavior. However, this explanation 
would appear to be inconsistent with the experiments of 



Babcock et al |107j discussed above, as well as with the 
evidence that hole doping near the grain boundaries by 
Ca substituting Y improves the conduction across grain 
boundaries. Klie et al argue that their calculations [T^T] 
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support the hole depletion scenario (hence positive po- 
tential at the grain boundary) and that this raises ques- 
tion about what potential is measured in the electron 
holography experiments. However, Koch'122' argues 
that "it seems to be possible to rule out that this obser- 
vation is a measurement artifact inherent to the method 
of electron holography" . 

These seemingly contradictory observations however 
are consistent with the physics predicted by the model 
under consideration here. We model the grain bound- 
ary by assuming a smaller hopping amplitude between 
sites at opposite sides of the grain boundary. A large 
angle grain boundary would presumably have a more 
reduced hopping amplitude compared to a small angle 
grain boundary. Figure 14 shows the charge distribution 
resulting from our model assuming a 70% reduction in 
the hopping amplitude across the grain boundary (de- 
noted by a dashed line), presumably corresponding to 
a fairly large angle grain boundary with significant in- 
crease in the resistance. The hole density in the vicinity 
of the grain boundary is significantly reduced, driven by 
lowering of kinetic energy of the system. As a conse- 
quence, the negative charge density becomes larger near 
the grain boundary, and this would give rise to a neg- 
ative electric potential at the grain boundary, consistent 
with the electron holography experiments, and a depleted 
hole concentration around the grain boundary, consistent 
with the EELS experiments [107] ■ 

As a function of temperature, Schneider et al |109j and 
Mennema et al |123| find that increasing the temperature 
leads to a rapid decrease of the grain boundary resis- 
tance. We suggests that this results from an increase in 
the hole concentration near the grain boundary as the 
temperature is increased. Figure 15 shows that as the 
temperature is increased the hole density near the grain 
boundary increases substantially in the dynamic Hub- 
bard model, and as a consequence the conductivity of the 
region would increases. It would be interesting to mea- 
sure directly the dependence of hole depletion on tem- 
perature by EELS experiments. This has not yet been 
done to our knowledge. 

Recent measurement of properties of grain boundaries 
in Fe-As superconductors [123] indicate that the behav- 
ior is similar to the high Tc cuprates. This is not sur- 
prising if the physics of both classes of materials, includ- 
ing their superconductivity [T^S], is described by dynamic 
Hubbard models. 

We expect that negative charge expulsion from the 
interior of superconductors should also lead to nega- 
tive charging of vortex cores in the mixed state. Such 
effect [126j and associated antiferromagnetic order in the 
vortex cores 1 127] has been detected experimentally in 
YBa2Cu307. 
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FIG. 14: Effect of a grain boundary, indicated by the dashed 
line for the conventional and dynamic Hubbard models with 
th ~ 0.1, U = 2. We assume that the hopping amplitude is 
reduced by a factor 0.3 for sites on opposite sides of the grain 
boundary. The hole occupation is depleted in the vicinity 
of the grain boundary in both cases, however the effect is 
much larger and extends over a wider range for the dynamic 
Hubbard model than for the conventional one. Temperature 
is T = 0.02. 
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FIG. 15: Same as Fig. 14 with temperature T = 0.1 



IX. SHAPE EFFECTS 

It is interesting to consider the effect of the shape of 
the sample on the charge expulsion profile in the dynamic 
Hubbard model. Consider an ellipsoidal shape as shown 
in Figure 16. The sites near the surface at the regions 
of higher curvature, i.e. top and bottom, have somewhat 
smaller hole concentration than at the regions of lower 
curvature at the lateral surfaces. This is easy to under- 
stand: the sites near the surface in the regions of high 
curvature have slightly lower coordination than those in 
the regions of low curvature, hence the holes do not ben- 
efit so much from kinetic energy lowering and prefer to 
stay away from those regions. Thinking in terms of elec- 
trons instead of holes, it means the body expels more 
electrons to the top and bottom than to the sides. This 
should give rise to a higher electric potential near the 
sides than at the top and bottom, and a quadrupolar 
electric field with field lines starting at the lateral sides 
and ending at the top and bottom. This is precisely the 
type of electric field found in the superconducting state 
by solving the alternative London equations proposed to 
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FIG. 16: Effect of shape on the charge expulsion profile. 
th — 0.1, U — 2. In the ellipsoidal shape shown, for the 
conventional Hubbard model the charge occupation near the 
surface is similar near the top surface and the lateral surfaces. 
For the dynamic Hubbard model the hole concentration is 
somewhat higher near the lateral surfaces than near the top 
and bottom surfaces. For the origin of this effect, see text. 
Temperature is T = 0.02. 



describe the electrodynamics of superconductors within 
this theory [m [HE]. 

A qualitative way to understand this charge distribu- 
tion in the superconducting state is the following: the 
electrons near the lateral surfaces can move faster than 
those at the top and bottom in the region of high cur- 
vature, just as racing cars. Hence they will have higher 
kinetic energy and consequently lower potential energy 
than the electrons near the top and bottom surfaces, so 
as to keep the same sum of kinetic and potential energies. 
Electrons near the lateral surfaces having lower potential 
energy means that the electric potential is higher near the 
lateral surfaces than near the top and bottom surfaces, 
resulting in electric field lines starting from the side and 
ending at the top and bottom just as found from the 
analysis of the hole distribution in the dynamic Hubbard 
model discussed in the above paragraph. 

More generally, using the same arguments we expect 
that for other body shapes the electric potential will be 
higher in the regions of lower surface curvature and lower 
in the regions of higher curvature in the dynamic Hub- 
bard model and in superconducting bodies. An example 
for a body shape resulting from combining a prolate and 
an oblate ellipsoid is given in Ref. |129) . 



X. CHARGE EXPULSION IN THE 
SUPERCONDUCTING STATE 

As seen in the previous sections, the dynamic Hubbard 
model has a tendency to expel negative charge from its 
interior to the surface driven by lowering of kinetic en- 
ergy. Starting with a charge neutral system in the normal 
state, where a uniform positive ionic charge distribution 
is compensated by an equal uniform negative electronic 
charge distribution, the negative charge expulsion would 
result in a net charge distribution as qualitatively shown 
in Fig. 2: a net excess positive charge in the interior and 
net excess negative charge near the surface. According to 
the numerical results in the previous sections (e.g. Fig. 
4) the positive charge in the interior predicted by the 
dynamic Hubbard model Hamiltonian is approximately 
uniform, just as that predicted from the electrodynamic 
equations in the superconductor [25] . 

This would result in the presence of an electric field 
in the interior of the system, that increases linearly in 
going from the center towards the surface. However, this 
cannot happen in a real normal metal since a metal in 
the absence of current flow cannot have a macroscopic 
electric field in the interior. Therefore, we conclude that 
longer range Coulomb interactions, omitted in the dy- 
namic Hubbard model, prevent this from actually taking 
place in a real material. In other words, potential en- 
ergy triumphs over kinetic energy in the normal state, 
and a macroscopic metal will remain charge neutral in 
the interior, despite the tendency to develop this macro- 
scopic charge inhomogeneity if dynamic Hubbard model 
physics is dominant. At most, the system will develop 
local charge inhomogeneity that will be screened within 
a Thomas Fermi length, that can be several A in sys- 
tems like underdoped high Tc cuprates where the carrier 
density is very low. 

However, the situation can change if the system enters 
the superconducting state at low temperatures. There 
is no a-priori reason why a superconductor cannot have 
an electric field in its interior [25]. A superconductor 
is a macroscopic quantum system, and quantum sys- 
tems in their ground state minimize the sum of poten- 
tial and kinetic energies. That should not in general re- 
sult in a uniform charge distribution that minimizes po- 
tential energy only. The electrodynamic equations that 
we have proposed for superconductors [2S] predict that 
the superconductor has rigidity in the charge degrees of 
freedom [25l llOOj and will not screen an interior electric 
field as a normal metal would. 

To compute the charge distribution in the super- 
conducting state we solve numerically the Bogoliubov 
de Gennes (BdG) equations for the dynamic Hubbard 
model, for systems with the same geometry as discussed 
in the previous sections. For the correlated hopping 
model {At2 = 0) the equations are given in Ref. |103| . 
and are simply extended for the case At2 7^ 0. There 
are two gap parameters. An and Aij corresponding to 
on-site and nearest-neighbor pairing amplitudes [TOS]. 
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FIG. 17: Comparison of occupations in the normal and su- 
perconducting states. Radius of cylinder is i? = 9, average 
occupation per site is n = 0.126. th = 0.1, U = 2, At = 0.2, 
At2 = 0.4, ksT — 0.02. The average gap parameters in 
the superconducting state are A;; = 0.0044, Aij — —0.025. 
We also show the occupations for the conventional Hubbard 
model. At = = At2 — 0. The charge expulsion is largest in 
the normal state of the dynamic Hubbard model. 



FIG. 19: Dynamic Hubbard model, parameters as in Fig. 17. 
Comparison of occupations right above Tc, ksT — 0.05, and 
at ksT = 0.02, with the system in the normal (full line) and 
superconducting (dashed line) state. Note that when the sys- 
tem is allowed to go superconducting, the occupations below 
Tc essentially don't change, while if kept in the normal state 
the charge expulsion increases as the temperature is lowered. 
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FIG. 18: Real space hole occupation for the cases of Fig. 
17 with At 7^ 0. Note that the interior hole occupation is 
slightly larger in the normal state, and the occupation near 
the surface slightly smaller. 



We have tested our computer program by solving the 
BdG equations numerically on a square lattice with peri- 
odic boundary conditions and comparing with the stan- 
dard BCS solution. For the cylindrical geometry with 
open boundary conditions considered here, the numerical 
solution obtained for the gap parameters deep in the in- 
terior are close to the gap parameters found in the square 
lattice with periodic boundary condition using both the 
BdG equations and the standard BCS equations applica- 
ble to translationally invariant systems. The gap param- 
eters go to zero as the surface is approached. 

Initially we had hoped[21] that comparison of the oc- 
cupations in the dynamic Hubbard model in the normal 
and superconducting states would yield clear evidence 
that the system expels negative charge from the interior 
to the surface as it goes superconducting, as expected on 



physical grounds [24 l 199 ] and predicted by the electrody- 
namic equations [3S]. This is not what happens, as shown 
in Figs. 17 and 18. Instead, the charge distribution be- 
comes more uniform in the superconducting compared 
to the normal state at the same temperature. In fact, 
it appears that as the temperature is lowered and the 
system enters into the superconducting state the charge 
expulsion that increases in the dynamic Hubbard model 
as the temperature is lowered in the normal state stops 
changing and stays essentially the same as what it is at 
Tc when the system is cooled below Tc- This is shown in 
Figure 19. 

In summary, from the numerical results presented here 
it appears that the BCS/BdG solution of the dynamic 
Hubbard model does not reflect the charge expulsion pre- 
dicted by the electrodynamic equations as the system 
enters the superconducting state ^25^. On the other hand 
this is perhaps not too surprising. The charge expulsion 
predicted by the electrodynamic equations is of the order 
of 1 extra electron every 10^ sites near the surface[13Q], 
which certainly would not be noticeable in systems of the 
size considered here. 



XI. SUPERCONDUCTIVITY AND CHARGE 
IMBALANCE 

Within the BCS formalism, the total electronic charge 
per site is given by 

Qtot = I Y.^ulf{Ek) +vl{l- fiEk))] (26) 

k 

in units of the charge of one carrier, e or — e depending 
on whether one is using electron or hole representation. 
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Eq. (26) can be written as |131j 

Qtot = Qc + Q* 



with 



the charge of the condensate, and 



(27) 



(28a) 



(28b) 



the average charge of the quasiparticles. The coherence 
factors are given by the usual form 



2 , efe-ptx 



(29a) 




FIG. 20: Gap function and quasiparticle energy Ek as a 
function of band energy in hole representation. 



which is given approximately by (again ignoring band 
edge effects) 



(1 



(29b) 



In a conventional BCS superconductor Q* = in equilib- 
rium since quasiparticles are charge neutral on average, 
half electron, half hole. A non-zero Q*, termed "charge 
imbalance" or "branch imbalance" , can be generated in 
a non-equilibrium situation in the presence of current 
flow[ 132t il33j and/or a temperature gradient [l341 1135j . 

Instead, in the dynamic Hubbard model (or the corre- 
lated hopping model) a branch imbalance exists even in 
equilibrium [60j . The gap function has a slope of universal 
sign [55] 



AfeEE A(efe) = A,„(-efe/(I?/2)-Kc) 



(30) 



with D the bandwidth and A„i > and c obtained from 
solution of the BCS equations [SS]. The minimum gap 



is Aq = A(/.t)/a, with a 
quasiparticle energy is given by 



{1\^I(PI2)Y and the 



(31) 



The minimum gap Aq is attained not at = /i but at 
^k = 1^ + V , with 



= A„,(T)Ao(r)/(ai?/2) >0. 



(32) 



Both Aq and A™ go to zero at Tc as \JTc — T so v goes 
to zero linearly as T approaches from below. The gap 
function and quasiparticle excitation spectrum are shown 
schematically in Figure 20 in hole representation. In 
equilibrium, quasiparticles are symmetrically distributed 
around the minimum located ai e'^ = fj, + v and as a con- 
sequence Q* > 0, quasiparticles are positively charged on 
average. If we ignore band edge effects we have simply 



^ N ^ Ek 



(33) 



-AT) 



-/3Ao 



Da (/3An)V2' 



(34) 



Q* is suppressed at low temperatures due to the expo- 
nential factor, peaks somewhat below Tc and goes to zero 
at Tc- Numerical examples are shown in ref. |60) . How- 
ever when the finite bandwidth is taken into account Q* 
remains positive at and above. 

Figure 21 shows the distribution as a function of the 
distance to the center r of the quasiparticle charge Q* at 
site I, given by 



Qt 



-Ti 



^)I{En 



(35) 



obtained from solution of the BdG equations, for a dy- 
namic Hubbard model and for an attractive Hubbard 
model. In Eq. (35), Uni and Vni are the amplitudes 
of the n-th eigenvector at site i obtained from diagonal- 
ization of the BdG Hamiltonian |103j . En is the energy 
for state n and / is the Fermi function. In the attrac- 
tive Hubbard model (Fig. 21 lower panel) particle-hole 
symmetry is broken only because the band is not half- 
full, but the interaction is particle-hole symmetric. As a 
consequence, the quasiparticle charge oscillates between 
positive and negative values. Instead, as seen in Fig. 21 
(upper panel) in the dynamic Hubbard model quasipar- 
ticles are predominantly positively charged, as expected 
due to the shift in the chemical potential by u displayed 
in Fig. 20. 

Figure 22 shows the real space distribution of the 
quasiparticle charge in the dynamic Hubbard model 
(right panel). The total site occupation for this case is 
shown on the left panel. It can be seen from Figs. 20 
and 21 that the positive quasiparticle charge is located 
mostly near the surface of the system. This is relevant 
to the discussion in the next section. 
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FIG. 21: Quasiparticle charge (Eq. (35)) as function of dis- 
tance to the center for a dynamic Hubbard model with param- 
eters as in Fig. 17 and for an attractive Hubbard model with 
th = 0.1, U = —0.4. Average occupation is n = 0.126, tem- 
perature is fcsT — 0.01. Note that in the dynamic Hubbard 
model the quasiparticle charge is predominantly positive. 
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FIG. 22: Real space hole occupation (left panel) and quasi- 
particle charge (right panel) in the superconducting state of 
the dynamic Hubbard model. Parameters are the same as for 
Figure 21 (upper panel). Note that the quasiparticle charge 
is predominantly near the surface of the sample. 



XII. TWO-FLUID MODEL AND INTERIOR 
ELECTRIC FIELD 



In this section we discuss to what extent the dynamic 
Hubbard model reflects the physics shown in Figure 2 in 
the superconducting state, how it depends on tempera- 



ture, and how this physics could be detected experimen- 
tally. 

We assume a two-fluid model, with the total carrier 
concentration independent of temperature. We have then 



Us = ns{T) + n„(T) 



(36) 



with Us the total carrier (hole) concentration and ns{T) 
and nn{T) the superfluid and normal components at 
< T < Tc. Within the two-fluid interpretation of BCS 
theory they are given in terms of the London penetration 
depth as 



ns{T) 



'Ai(T) 



r) 



(37a) 



(37b) 



with Xl (Al(T')) the London penetration depths at zero 
(finite) temperature. Ignoring finite bandwidth effects 
the quasiparticle density is then [136 

which yields at low temperatures 

n„(r) = V2^(/3Ao)i/2e-^^°n,. (39) 

On the other hand, the average quasiparticle charge per 
site is given by Eq. (34). Combining with Eq. (39), 



Q*(r) IfcsT 



knTAr. 



,(T) aD/2Ao Al + {D/2)^ ^ {D/2y 



(40) 



It can be seen that the quasiparticle charge is a small 
fraction of the quasiparticle density. For example, for 
the parameters of Fig. 21 we have 



D = 2z{th + nAt + ^Ata) 1-01, 



(41) 



A,j = 0.00253, A™ = zAy = 0.10, hence Q*/n„ = 
0.004. 

Assuming the system as a whole is charge neutral, the 
negative charge of the electrons in the band exactly com- 
pensates the positive charge of the ions, which is uni- 
formly distributed in space (except for variations on the 
scale of A). At temperatures below T^, the quasiparti- 
cles have a net positive charge, hence as a consequence 
the condensate has a total negative charge greater than 
the total positive charge of the ions. The condensate is 
highly mobile, and just like in a normal metal any excess 
negative charge will move to the surface [137] it is natural 
to expect that negative charge from the condensate will 
move to the surface. 

Furthermore, we have seen in the previous section that 
the positive quasiparticle charge is located predominantly 
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near the surface in the superconducting state (Fig. 22 
right panel). This can be interpreted as reflecting the 
fact that the superfluid has higher negative density near 
the surface, and the positive normal fluid consequently 
develops higher density near the surface to screen the 
superfluid charge. In addition, as already seen in the 
normal state of the dynamic Hubbard model, the total 
hole concentration is smaller near the surface which im- 
plies extra negative charge near the surface. Thus we ar- 
gue that the dynamic Hubbard model provides support 
to the prediction of the electrodynamic equations of the 
theory[531 [25] that the superconductor expels superfluid 
negative charge from the interior to the surface. 

Whether or not a macroscopic electric field will exist 
in the interior of the superconductor depends on whether 
there are enough quasiparticles to screen the electric field 
created by the negative charge expulsion of the conden- 
sate. In the ground state (no quasiparticles) the theory 
predicts that the net positive charge density in the inte- 
rior is |130| 



Pa 



Ps + Pi, 



(42) 



with R the radius of the cylinder and — h/{2m^c) — 
O.OOI93A the quantum electron radius, and there is a 
negative charge density 



P- 



R 
2A 



-PO^ Ps+ Pi. 



(43) 



within a London penetration depth of the surface, as 
shown schematically in Fig. 2. The charge densities at 
zero temperature are shown schematically in the left pan- 
els of Figure 23. ps denotes the superfluid charge density. 

At finite temperatures, there is a positive quasiparti- 
cle charge density excited, p„ = |e|(5*, and a crossover 
temperature Tcr can be defined. For temperatures lower 
than Ten the average quasiparticle charge density excited 
is smaller than pq . In the middle top panel of Fig. 23 we 
assume pn is uniformly distributed, and in the right top 
panel we assume all p„ has moved to within the London 
penetration depth of the surface. Even so, it is unable 
to screen the internal electric field, since a positive net 
charge density p ^ po — Pn > remains in the interior 
and a negative net charge density p < remains near the 
surface, as shown in the top right panel of Fig. 23. At 
the crossover temperature T^r the quasiparticle charge 
density excited reaches the value po, and by migrating 
to the region within a London penetration depth of the 
surface (lower right panel of Fig. 23) it can completely 
screen both the interior positive charge and the negative 
charge in the surface layer, so that the electric field ev- 
erywhere gets cancelled. This will also be the case at any 
temperature T > Tcr- 

The value of the crossover temperature can be obtained 
from the equation 



Q* 



Q* 



rir, 



R 



(44) 
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FIG. 23: Charge distribution in the superconductor 
(schematic). In the ground state (left panels) there is an 
excess positive charge density po in the interior and an excess 
negative charge density p_ near the surface. At low temper- 
atures T < Tcr (upper panels) the charge density is partially 
screened by excited quasiparticles but a net positive charge 
density p > remains in the interior and a net negative charge 
density p < near the surface (dot-dashed lines in upper right 
panel). At temperature T — Tcr and above the net charge 
density is zero both in the interior and near the surface. 



with Q* /un given by Eq. (40) and n„ given by Eq. (39), 
hence 



ksTcr^m , 27rAo 
(-0/2)2 ^ksTcr 



)l/2g-Ao/(fcBT,,) 



For example, assuming the usual relation 2Ao/fcsrc — 
3.53, for the case under consideration with (5*/n„ — 
0.004 and assuming R ^ 500nm yields T^,. = 0.16Tc. 
For temperatures lower than T^, a nonzero electric field 
is predicted to exist in the interior of the superconductor. 

In the presence of a non-zero internal electric field, 
superconductors of non-spherical shape should also de- 
velop electric fields extending to the region exterior to 
the body, of magnitude and direction determined by the 
shape of the body and the electrodynamic equations of 
the superconductor [1^51 1129 j . These electric fields should 
be experimentally detectable in the neighborhood of su- 
perconductors at temperature lower than T^r- In ad- 
dition, the internal electric field should be directly de- 
tectable in electron holography experiments, as discussed 
in a recent paper [138] . The magnitude of these predicted 
electric fields is of order of Hd, the lower critical mag- 
netic field, in the interior of the superconductor [1^0] . and 
an appreciable fraction of it in the region outside the su- 
perconductor near the surface, depending on the shape of 
the bodv |12gl 1129] . No external electric field is expected 
outside a planar surface or a spherical body. 



17 



XIII. THE MEISSNER EFFECT, THE LONDON 
MOMENT AND THE GYROMAGNETIC EFFECT 

The fact that in superconductors the superfluid carries 
negative charge is estabhshed experimentally from exper- 
iments that measure the London moment |139[ 1140] : a ro- 
tating superconductor develops a magnetic field that is 
always parallel^ never antiparallel, to the direction of the 
mechanical angular momentum of the body [55]. 

We have seen that the dynamic Hubbard model has a 
tendency to expel negative charge, that in a real system is 
inhibited in the normal state because of the effect of long- 
range Coulomb repulsion but may take place when the 
system becomes superconducting. The considerations in 
the previous section suggest that as a system becomes 
superconducting the electrons that condense into the su- 
perfluid state are partially expelled towards the surface, 
and at the same time normal electrons flow inward to 
compensate for the charge imbalance, as indicated by the 
fact that the positive quasiparticle charge moves outward 
in the superconducting state as seen in the last section. 

Consider now these processes in the presence of an 
external magnetic field in the z direction, as shown in 
Fig. 24. The outflowing superfluid electrons will be de- 
flected counterclockwise by the Lorentz force exerted by 
the magnetic field, building up a Meissner current flow- 
ing clockwise near the surface that suppresses the applied 
field in the interior. At the same time, the inflowing 
normal electrons are defiected clockwise by the magnetic 
field. Because these electrons undergo scattering from 
the ions, they will transmit their momentum to the ions 
and the body as a whole will start rotating in a clock- 
wise direction. And because these electrons are slowed 
down and ultimately stopped by the collisions with the 
ions they will not reinforce the applied magnetic field. 
The end result is a superfluid current near the surface 
flowing in clockwise direction (i.e. superfluid electrons 
flowing in counterclockwise direction) that suppresses the 
interior magnetic field, and a slow body rotation in the 
clockwise direction that exactly cancels the mechanical 
angular momentum carried by the superfluid electrons in 
the Meissner current, as required by angular momentum 
conservation. 

The magnetic held generated by rotating supercon- 
ductors (London moment) [1391 I140j can be similarly 
explained '141' by the fact that in a rotating normal metal 
that is cooled into the superconducting state the expelled 
superfluid electrons, that are moving at the same angu- 
lar velocity as the body, will have a smaller tangential 
velocity than the body when they reach the surface, giv- 
ing rise to a net current and resulting magnetic moment 
in direction parallel, never antiparallel, to the angular 
momentum of the rotating body. 

More quantitatively, the outflow occurs because super- 
fluid electrons expand their orbits from microscopic ra- 
dius kp^ to mesoscopic radius 2Al|142I, in the process 



normal electron 




FIG. 24: The outflowing superfluid electrons are deflected in 
the counterclockwise direction by the applied magnetic fleld, 
generating a clockwise Meissner current (iMeissner) near the 
surface that suppresses the magnetic field in the interior. The 
inflowing normal electrons collide with the ions and impart 
the entire body with angular momentum antiparallel to the 
applied magnetic field, that is equal and opposite to the me- 
chanical angular momentum of the electrons in the Meissner 
current . 

acquiring an azimuthal velocitv |143| 



which is the speed of the superfluid electrons in the Meiss- 
ner current [136] . The total mechanical angular momen- 
tum acquired by these electrons in a cylinder of radius R 
and height h is 

Lei = iTR^hns{meV^2\L) (47) 

which coincides with the total electronic angular momen- 
tum of the Meissner current flowing within a London pen- 
etration depth of the surface 

LMeissner = 27ri?AL /in^ (TOgW^i?) . (48) 

The inflowing normal electrons transmit the same an- 
gular momentum to the body as a whole by colli- 
sions with the ions, as required by angular momentum 
conservation [144] . As a consequence, the body starts 
rotating with angular velocity determined by the con- 
dition of angular momentum conservation, as measured 
experimentally [HSHIIZ] (gyromagnetic effect). 

These processes provide an intuitive explanation for 
the dynamics of the Meissner effect [T^ . the genera- 
tion of the London moment, the gyromagnetic effect and 
the puzzle of angular momentum conservation [1441 1148j 
in superconductors, provided the superconductor is de- 
scribed by a dynamic Hubbard model that gives rise 
to negative charge expulsion. In superconductors not 
described by dynamic Hubbard models but by conven- 
tional BCS-electron-phonon theorv [T36] no charge ex- 
pulsion takes place so these considerations don't apply. 
In those superconductors the dynamical origins of the 
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Meissner effect and the London moment and the expla- 
nation of the angular momentum puzzle remain to be 
elucidated. 



XIV. DISCUSSION 

Both the conventional Hubbard model and the dy- 
namic Hubbard model are simplified descriptions of real 
materials, and whether or not they contain the physics 
of interest for particular real materials is in principle an 
open question. In this paper we have argued that the 
new physics that the dynamic Hubbard model incorpo- 
rates beyond what is contained in the conventional Hub- 
bard model has vast implications for the understanding 
of many materials of current interest. 

The new physics of the dynamic Hubbard model is 
that it allows the electronic orbital to expand when it is 
doubly occupied. This expansion has associated with it 
outward motion of negative charge as well as lowering of 
the electron's kinetic energy at the atomic level, and it 
is of course not "electron-hole symmetric" . This physics 
is not included in the conventional Hubbard model that 
assumes that the electronic orbital is infinitely rigid. The 
key question is not whether this physics exists in real 
atoms, of course it does 3 , the key questions are how 
large the effect is, what are its consequences, and under 
which conditions and for which materials is it or is it not 
relevant for various properties of the materials. 

The mathematical treatment of the model is not sim- 
ple, and from the outset we restricted ourselves in this 
paper to the antiadiabatic limit, i.e. assuming that the 
energy scale associated with the orbit expansion (ujq in 
Eq. (4)) is sufficiently large than it can be assumed 
infinite. This brings about the simplification that the 
high energy degrees of freedom can be eliminated and 
the Hamiltonian becomes equivalent to the low energy 
effective Hamiltonian Eq. (7), a Hubbard model with 
correlated hoppings, linear and nonlinear terms At and 
At2. This low energy effective Hamiltonian, together 
with the quasiparticle weight renormalization described 
by Eq. (6), describes many properties that we believe 
are relevant to real materials and are not described by 
the conventional Hubbard model. In particular it gives 
rise to hole superconductivity [5S], driven by lowering of 
kinetic energy of the carriers [5^ 1 149j . with many charac- 
teristic features that resemble properties of the cuprates. 
In other work we have also examined the effect of the high 
energy degrees of freedom in describing spectral weight 
transfer from high to low energies ( "undressing" [86 ) as 
the number of holes increases and as the system enters 
the superconducting state, as well as the effect of finite 
ujQ in further promoting pairing in this modeljSSj. 

In this paper we have shown with specific quantitative 
examples that the dynamic Hubbard model has a strong 
tendency to expel negative charge from the interior of 
the system to the surface, driven by lowering of kinetic 
energy of the charge carriers. We believe it is truly no- 



table that this property of the model mimics the physics 
of the single atom that motivated the formulation of the 
model, even though the kinetic energy lowering and nega- 
tive charge cloud expansion of the atomic electron is not 
explicitly included in the site Hamiltonian Eq. (1). It 
is furthermore notable that the orbital expansion in an 
ion is larger when the ion is negatively charged, which 
in the model translates to a larger value of the coupling 
constant g (smaller Eq. (6b)), larger At/th (Eq. (10)) 
and consequently larger tendency for the model to expel 
negative charge. We argue that the fact that the lattice 
Hamiltonian Eq. (9) describes at a macroscopic length 
scale the same physics at the atomic scale that motivated 
the Hamiltonian is a strong indication that the Hamilto- 
nian is relevant for the description of reality. 

The tendency of the dynamic Hubbard model to expel 
negative charge and the tendency to pairing of holes and 
superconductivity of course go hand in hand: they both 
originate in the fact that the kinetic energy of a hole 
is lowered when another hole is nearby. Increasing the 
chance of having another hole nearby can be achieved 
by negative charge expulsion, thus increasing the overall 
hole density, and by pairing, thus increasing the local 
hole density. 

The effects predicted by this Hamiltonian are largest 
when the coupling constant g is large, or equivalently 
when the overlap matrix element S is small, which cor- 
responds to a "soft orbital" that would exist for neg- 
atively charged anions, and the effects are also largest 
when the band is almost full with negative electrons 
(strong coupling regime) |150j . Thus, not surprisingly, 
more negative charge at the ion or/and in the band 
yield larger tendency to negative charge expulsion for 
the entire system. We believe that the Hamiltonian is 
relevant to describe the physics of many superconduc- 
tors of current interest including high cuprates, Fe 
pnictides, Fe chalcogenides, MgB2 and BiS-)-haseA \lb\\ 
superconductors. All these materials have negatively 
charged ions {O , As ,5 ,Se ,-B^) with soft 
orbitals, and for most of them, including "electron- 
doped" cuprates |152| . there is experimental evidence for 
dominant hole transport in the normal state. We sug- 
gest that the orbital expansion and contraction of these 
negative ions depending on their electronic occupation is 
responsible for many interesting properties of these ma- 
terials including their superconductivity, and is described 
by the dynamic Hubbard model. 

In addition, we have seen in this paper that the model 
leads to charge inhomogeneity driven by lowering of ki- 
netic energy, and in extreme cases to phase separation, 
and that it leads to negatively charged grain boundaries 
and depletion of hole carriers in the vicinity of grain 
boundaries, properties experimentally observed in many 
of these materials but not understood using conventional 
models such as "band bending" . 

We have also examined here the question whether the 
interiorpS] and exterior [128| electric fields predicted to 
exist in the ground state of superconductors within this 
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theory would exist also at finite temperatures and con- 
cluded that they should exist and hence be experimen- 
tally detectable up to a crossover temperature T^r, cal- 
culated to be 0.1 6Tc in one example. 

Finally, we have proposed that the negative charge ex- 
pulsion predicted by dynamic Hubbard models is relevant 
to the understanding of the Meissner effect, the London 
moment and the gyromagnetic effect exhibited by super- 
conductors that are not described by conventional BCS- 
electron phonon theory. 

In summary, we hope the reader will appreciate the re- 
markable qualitative similarity of Figs. 1 and 2, depicting 
an atom and a superconductor within this model. Super- 
conductors have been called "giant atoms" in the early 
days of superconductivity for many other reasons [1531 - 
I155| . The essential property of the atom, that it is not 
electron-hole symmetric because the negative electron is 
much lighter than the positive nucleus, manifests itself 
identically in the atom described by the dynamic Hub- 
bard model and in the state of a macroscopic supercon- 
ducting body described by the model, and is missed in the 
world described by particle-hole symmetric conventional 



Hubbard models both at the atomic level and at the level 
of the macroscopic superconductor. The superconductor 
is indeed a "giant atom" within our description, with the 
highly mobile light negative superfluid reflecting the nim- 
ble atomic electron, and the heavy positive quasiparticles 
reflecting the positive nucleus. 

Much of the physics of dynamic Hubbard models for 
finite wq remains to be understood. In fact, the model 
itself may require substantial modification to account for 
different values of ujq for different electronic occupations: 
the excitation spectrum of H is certainly very different 
from that of H~ . In connection with this and going be- 
yond the antiadiabatic limit where only diagonal transi- 
tions of the auxiliary boson field are taken into account 
as in this paper, it is possible that vertical transitions 
may play a key role in describing the superconducting 
state [156] . It is also an open question to what extent 
dynamic Hubbard models can describe the mysterious 
"pseudogap state" of underdoped high Tc cuprate mate- 
rials. These and other questions will be the subject of 
future work. 
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